

\section{有限元基础}

\def\s1eq{有限元基础：}

\subsection{边值问题}
典型的边值问题可用区域$\Omega$内的控制方程和包围区域$\Omega$的边界$\varGamma$上的边界条件来定义。
控制方程可表示为:
\begin{equation} \label{s1eq 控制方程}
    \mathcal{L} \varPhi =f
\end{equation}
$\mathcal{L}$是算子，$f$是激励或强加函数，$\varPhi$是代求未知量。

\subsection{里兹变分法}
里兹法是一种变分方法，其中边值问题可用其对应的泛函极值表示，泛函极小值即为给定边界条件下
的控制方程的解。通过求泛函相对于其变量的极小值，可得到控制方程的解。
\begin{definition}
    (内积)
    \begin{equation}
        \left< \varPhi ,\varPsi \right> =\int_{\varOmega}{\varPhi \varPsi ^*d\varOmega}
    \end{equation}
    星号*表示复共轭。
\end{definition}
\begin{definition}
    (自伴算子)
    若算子满足：
    \begin{equation}
        \left< \mathcal{L} \varPhi ,\varPsi \right> =\left< \varPhi ,\mathcal{L} \varPsi \right> 
        \end{equation}
        则称算子$\mathcal{L}$是自伴的。
    \end{definition}
\begin{definition}
    (正定算子)
    若算子满足：
    \begin{equation}
        \left< \mathcal{L} \varPhi ,\varPhi \right> = \left\{ \begin{array}{c}
            \begin{matrix}
            >0&		\varPhi \ne 0\\
        \end{matrix}\\
            \begin{matrix}
            =0&		\varPhi =0\\
        \end{matrix}\\
        \end{array} \right. 
    \end{equation}
则算子$\mathcal{L}$是正定的。
\end{definition}
若算子$\mathcal{L}$既是自伴又是正定，则控制方程的解即为
下面泛函对$\widetilde{\varPhi}$的极小值：
\begin{equation} \label{\s1eq 控制方程对应的泛函}
    F\left( \widetilde{\varPhi } \right) =\frac{1}{2}\left< \mathcal{L} \widetilde{\varPhi },\widetilde{\varPhi } \right> -\frac{1}{2}\left< \widetilde{\varPhi },f \right> -\frac{1}{2}\left< f,\widetilde{\varPhi } \right> 
\end{equation}
式中，$\widetilde{\varPhi }$为试探解或称为近似解。

通过下面步骤求上述泛函对$\widetilde{\varPhi}$的极小值，以下讨论中，
假设问题均为实数值,则可假设$\widetilde{\varPhi}$的近似展开为：
\begin{equation} \label{\s1eq 试探解的近似展开}
    \widetilde{\varPhi }=\sum_{j=1}^N{c_jv_j}=\left\{ c \right\} ^T\left\{ v \right\} =\left\{ v \right\} ^T\left\{ c \right\} 
\end{equation}
式中，$v_j$是定义在全域上的展开函数，$c_j$是待定展开系数。

将式(\ref{\s1eq 试探解的近似展开})待入式(\ref{\s1eq 控制方程对应的泛函})中，可以得到：
\begin{equation}
F=\frac{1}{2}\left\{ c \right\} ^T\int_{\Omega}{\left\{ v \right\} \mathcal{L} \left\{ v \right\} ^T}d\Omega \left\{ c \right\} -\left\{ c \right\} ^T\int_{\Omega}{\left\{ v \right\} fd}\Omega 
\end{equation}

令$F$对$c_i$的偏导为零，有：
\begin{equation}
    \begin{aligned}
        \frac{\partial F}{\partial c_i} & =\frac{1}{2}\int_{\Omega}{v_i\mathcal{L} \left\{ v \right\} ^T\mathrm{d}\Omega \left\{ c \right\}}+\frac{1}{2}\left\{ c \right\} ^T\int_{\Omega}{\left\{ v \right\} \mathcal{L} v_i\mathrm{d}\Omega}-\int_{\Omega}{v_if\mathrm{d}\Omega}
        \\
        &=\frac{1}{2}\sum_{j=1}^N{c_j}\int_{\Omega}{\left( v_i\mathcal{L} v_j+v_j\mathcal{L} v_i \right) \mathrm{d}\Omega}-\int_{\Omega}{v_if\mathrm{d}\Omega}
        \\
        &=0\begin{matrix}
            &		&		i=1,2,3,\cdots ,N\\
        \end{matrix}       
    \end{aligned}
\end{equation}
从而得到线性代数方程组，为了更清晰表示其形式，通过变量替换，可以得到下面矩阵方程：
\begin{equation}
    \left[ S \right] \left\{ c \right\} =\left\{ b \right\} 
\end{equation}
\begin{note}
    这里[*]表示矩阵，\{*\}表示向量。
\end{note}
\noindent 通过求解该矩阵方程，可以得到控制方程(\ref{s1eq 控制方程})的近似解。

\noindent $[S]$的元素用下标表示法可以写成：
\begin{equation}
    S_{ij} = \frac{1}{2} \int_{\Omega}( v_{i}\mathcal{L}v_{j} + v_{j}\mathcal{L}v_{i} ) d\Omega
\end{equation}
$\{ b \} $的元素用下标表示法可以写成：
\begin{equation}
    b_{i} = \int_{\Omega} v_{i}f d\Omega
\end{equation}
由于算子$\mathcal{L}$是自伴的，因此$S_{ij}$是对称的，可以表示为：
\begin{equation}
    S_{ij} = \int_{\Omega}v_{i}\mathcal{L}v_{j} \mathrm{d}\Omega
\end{equation}

\subsection{伽辽金方法}
伽辽金方法属于残数加权方法类型，它通过对控制方程的残数求加权的方法来得到方程的解。若
设$\widetilde{\varPhi}$是控制方程(ref{\s1eq 控制方程})的近似解，由于解是近似的，因此控制方程
并不严格为零，即存在非零残数：
\begin{equation}
    r = \mathcal{L}\widetilde{\varPhi}-f \ne 0
\end{equation}
\begin{definition}
    (残数加权)
    为了使得近似解$\widetilde{\varPhi}$在$\Omega$内所有点上有最小值，定义所有点上残数加权的积分公式：
    \begin{equation} \label{\s1eq 残数}
        R_{i} = \int_{\Omega}\omega_{i} r \mathrm{d}\Omega
    \end{equation}
    式中，$R_{i}$表示所有点上残数加权的积分，$\omega_{i}$即为所选择的加权函数。
\end{definition}
通常加权函数有多种选择，常设加权函数与近似解展开式(\ref{\s1eq 试探解的近似展开})中所用的函数相同，即加权函数选择为：
\begin{equation}
    \omega _i=v_i
    \begin{matrix}
        &		&		i=1,2,3,\cdots ,N\\
    \end{matrix}
\end{equation}
则式(\ref{\s1eq 残数})变为：
\begin{equation}
    R_{i} = \int_{\Omega}( v_{i}\mathcal{L}\{v\}^{T}\{c\} - v_{i}f ) \mathrm{d}\Omega = 0
    \begin{matrix}
        &		&		i=1,2,3,\cdots ,N\\
    \end{matrix} 
\end{equation}
可以看到，伽辽金方法与里兹法得到了相同的方程组。

\subsection{有限元方法的基本步骤}
有限元方法是求边值问题的数值过程。该方法的原理是用许多子域来代表整个连续的区域。在子域中，
未知函数用带有未知系数的简单插值函数来表示。无限个自由度的原边值问题被转化成了有限个自由度
的问题，整个系统的解用有限数目的未知系数近似。然后，用里兹变分或伽辽金方法得到一组代数方程
（方程组）。最后，通过求解方程组得到原边值问题的解。综上所述，边值问题的有限元分析包含以下
步骤：\\
\indent（1）区域的离散或子域的划分；\\
\indent（2）插值函数的选择；\\
\indent（3）方程组的建立；\\
\indent（4）方程组的求解。\\

\subsubsection{区域离散}
区域$\Omega$被划分成许多小区域$\Omega^{e}(e=1,2,3,\dots,M)$表示，这些小区域称为子域或单元，
其中$M$表示子域总数。对于一维区域，子域是短直线段；二维区域，子域是三角形或矩形；三维区域
子域是四面体、三棱柱或矩形块。其中，线性线元、三角形单元以及四面体单元分别是基本的一维、二维
和三维单元。

多数有限元解中，问题是用与单元有关的结点上的未知函数$\varPhi$表示。线性单元有两个结点，对应
线段的两个端点；线性三角形单元有三个结点，对应三个顶点；线性四面体单元有四个结点，对应各四个
顶点。完整的描述一个结点需要知道它的坐标值、局部编码和全局编码。局部编码描述其在局部单元中的
编号，全局编码描述其在整个区域中的全局编号。
\begin{note}
    有限元法计算过程中，常得到带状矩阵，其带宽由一个单元中的两个结点的全局编码之差的最大值
    决定。通过适当的结点编码可以使得带宽最小，大大减少计算机的储存量和运算时间。

    区域离散作为一项预处理工作，可以独立于其它步骤。
\end{note}
\subsubsection{插值函数的选择}
插值函数用来近似表达一个单元中的未知解。通常，插值函数可以选择为一阶（线性）、二阶（二次）、或
高阶多项式。虽然更高的阶次有更高的精度，但通常会极大的增加计算量。在一般问题中，常采用一阶线性
插值，单元$e$结点对应的近似解可以表示为：
\begin{equation} \label{\s1eq 单元结点近似解}
    \widetilde\varPhi^{e} = \sum_{j=1}^{n} N_{j}^{e}\varPhi_{j}^{e}=\{N^{e}\}^T\{\varPhi^{e}\}
    =\{\varPhi^e\}^T\{N^e\}
\end{equation}
式中，$n$是单元中的结点数；$\varPhi_{j}^{e}$是单元$e$中结点$j$(此处为全局编号)的值；$N_{j}^{e}$是插值函数
或称为展开函数，其最高阶就是该单元的阶，插值函数的重要特征在于：它们只有在单元$e$内才不为零，
在单元$e$外均为零。

\subsubsection{方程组的建立}
1. 基于里兹法的求解公式

假设问题是实值的，则(\ref{\s1eq 控制方程对应的泛函})式给出的泛函$F$可以表示为：
\begin{equation}  \label{\s1eq 泛函展开式}
    F(\widetilde\varPhi) = \sum_{e=1}^{M}F^{e} (\widetilde\varPhi^{e})
\end{equation}
式中，$M$是全域的单元数，对于单元$e$，泛函$F^{e} (\widetilde\varPhi^{e})$有:
\begin{equation} \label{\s1eq 单元泛函表示}
    F^{e} (\widetilde\varPhi^{e}) = \frac{1}{2} \int_{\Omega^{e}} \widetilde{\varPhi}^{e} \mathcal{L} 
    \widetilde{\varPhi}^{e} \mathbf{d}\Omega - \int_{\Omega^{e}} f \widetilde{\varPhi}^{e} \mathbf{d}\Omega
\end{equation}
将式(\ref{\s1eq 单元结点近似解})带入式(\ref{\s1eq 单元泛函表示})，可以得到：
\begin{equation}
    F^{e} = \frac{1}{2} \{ \widetilde{\varPhi}^{e} \}^T \int_{\Omega^{e}} \{ N^{e} \}
    \mathcal{L} \{ N^{e} \}^T \mathbf{d}\Omega \{ \varPhi^{e} \} = \{ \varPhi^{e} \}^T
    \int_{\Omega^{e}} f \{ N^{e} \} \mathbf{d} \Omega
\end{equation}
将其写成矩阵形式：
\begin{equation} \label{\s1eq 单元泛函矩阵表示}
    F^{e} = \frac{1}{2} \{ \varPhi^{e} \} ^{T} [K^{e}] \{ \varPhi^{e} \} - \{ \varPhi^{e} \} ^{T}\{b^{e}\}
\end{equation}
式中，$[K^{e}]$是$n\times n$矩阵，$\{ b^{e} \}$是$n\times 1$列向量(这里$n$为单元$e$中包含的结点总数)：
\begin{equation}
    K_{ij}^{e} = \int_{\Omega^{e}} N_{i}^{e} \mathcal{L} N_{j}^{e} \mathbf{d} \Omega
\end{equation}
\begin{equation}
    b_{i}^{e} = \int_{\Omega^{e}} f N_{i}^{e} \mathbf{d}\Omega
\end{equation}
由于算子$\mathcal{L}$是自伴算子，因此，单元矩阵$[K^{e}]$是对称的。

将式(\ref{\s1eq 单元泛函矩阵表示})带入到式(\ref{\s1eq 泛函展开式})中，得到：
\begin{equation}
    F(\widetilde{\varPhi}) = \sum_{e=1}^{M} ( \frac{1}{2}\{ \varPhi^{e} \}^{T}
    [K^{e}] \{ \varPhi^{e} \} - \{ \varPhi^{e} \}^{T} \{b^{e}\} )
\end{equation}
上式存在求和符号，根据插值函数的特征：其只有在单元 e 内才不为零，在单元 e 外均为零。
将单元结点局部编码采用全局结点编码代替，上式可以写为：
\begin{equation}
    F = \frac{1}{2} \{ \varPhi \}^{T} [K] \{ \varPhi \} - \{ \varPhi \}^{T}\{b\}
\end{equation}
式中，$N$为全局结点总数，$[K]$是$N \times N$对称矩阵；$\{ \varPhi \}$是$N \times 1$
未知向量；$[b]$是$N \times 1$已知向量。为了求得未知量，令泛函$F$关于$\varPhi_{i}$的偏导为零：
\begin{equation}
    \frac{\partial F}{\partial \varPhi_{i}} = \frac{1}{2} \sum_{j=1}^{N} (K_{ij}+K_{ji})
    \varPhi_{j} - b_{i} = 0  
    \begin{matrix}
        & & & j = 1,2,3,\cdots,N
    \end{matrix}
\end{equation}
由于$[K]$是对称的，即$K_{ij} = K_{ji}$，因此，上式可以变为：
\begin{equation}
    \frac{\partial F}{\partial \varPhi_{i}} = \sum_{j=1}^{N} K_{ij}\varPhi_{j}-
    b_{i} = 0
    \begin{matrix}
        & & & j = 1,2,3,\cdots,N
    \end{matrix}
\end{equation}
写成矩阵形式有：
\begin{equation}
    [K]\{ \varPhi \} = \{ b \}
\end{equation}
上述推导过程是先对各个单元的泛函求和，然后再微分；也可以先对每个单元的泛函微分，然后再求和。

2. 基于伽辽金法的求解公式
第$e$各单元的残数加权为：
\begin{equation}
    R_{i}^{e} = \int_{\Omega_{e}} N_{i}^{e} (\mathcal{L} \widetilde{\varPhi}^{e}-f)
    \mathbf{d} \Omega
    \begin{matrix}
        & & & i = 1,2,3,\cdots,n
    \end{matrix}
\end{equation}
带入$\widetilde{\varPhi}^{e}$的近似解式(\ref{\s1eq 单元结点近似解})得到：
\begin{equation}
    R_{i}^{e} = \int_{\Omega^{e}} N_{i}^{e} \mathcal{L} \{N^{e}\}^{T} \mathbf{d}\Omega
    \{ \varPhi^{e} \} - \int_{\Omega^{e}} f N_{i}^{e} \mathbf{d} \Omega
    \begin{matrix}
        & & & i = 1,2,3,\cdots,n
    \end{matrix}
\end{equation}
写成矩阵形式为：
\begin{equation}
    \{ R^{e} \} = [K^{e}]\{ \varPhi^{e} \} - {b^{e}}
\end{equation}
利用局部编码和全局编码的关系，扩展上式得到：
\begin{equation}
    \{R\} = \sum_{e=1}^{M} \{\overline{R}^{e}\} = \sum_{e=1}^{M}( [\overline{K}^{e}]
    \{\overline{\varPhi}^{e}\} - {\overline{b}^{e}} )
\end{equation}
令上式等于零，可以得到方程组：
\begin{equation}
    \sum_{e=1}^{M} ( [\overline{K}^{e}]\{\overline{\varPhi}^{e}\} - {\overline{b}^{e}} ) = \{0\}
\end{equation}


\subsubsection{方程组的求解}
经过上述处理，最终得到的方程组的形式：
\begin{equation} \label{\s1eq 确定型方程}
    [K]\{ \varPhi \} = \{ b \}
\end{equation}
或者：
\begin{equation}  \label{\s1eq 本征型方程}
    [A]\{ \varPhi \} = \lambda[B] \{ \varPhi \}
\end{equation}
方程(\ref{\s1eq 确定型方程})是确定型的，它从非齐次线性微分方程或非齐次边界条件或从它们两者兼有的问题
中导出，在电磁学中，与存在源或激励的确定性问题有关。方程(\ref{\s1eq 本征型方程})是本征值型的，它从
齐次控制微分方程和齐次边界条件导出，在电磁学中，常与无源问题有关。

求解方程组得到$\{ \varPhi \}$，进一步计算就可以得到其它所需参数，并可以通过曲线、云图等形式
展示。这个步骤称为后处理，同区域离散一样可以同其它步骤分开。


